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Abstract. We show results of Monte Carlo simulations of an array of monodispersed 
magnetic monodomain particles, in a square lattice with dipolar interactions and perpendicular 
uniaxial anisotropy. We first show the equilibrium phase diagram of the system, 
which shows three phases, superparamagnetic, out-of-plane antiferromagnetic and in-plane 
antiferromagnetic with a reorientation transition between the last two. We then employ 
a recently introduced Time Quantified Monte Carlo method to study the relaxation of 
autocorrelations of the particles array for different ratios of dipolar to anisotropy energies. In 
the non-interacting case we show that relaxation is exponential in time with characteristic times 
obeying a classic result by Brown. When dipolar interactions are switched on, the relaxation 
is very well described by stretched exponential forms in the whole time window and in both 
the superparamagnetic and ordered phases. Relaxation times still obey a nearly Arrhenius 
behavior, with a single effective energy barrier that decreases as the dipolar interaction 
increases, a result that must be interpreted within the dynamics protocol. No signs of glassy 
behavior were found, in agreement with the absence of disorder in the model system. 
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1. Introduction 

The development of advanced experimental tools to fabricate magnetic nano-sctructures 
in a controlled way has led to an enormous growth of studies concerning the physical 
understanding of nanomagnetism, such as the role of different interactions and scaling 
behavior. Nano structured magnetic materials like magnetic particles and patterned magnetic 
alloys are very important systems in the experimental and theoretical research [1], beyond 
of its important technological applications like magnetic data storage and others fields 
(21 0. In general, magnetic media for data storage are composed of tiny but isolated 
magnetic nanocrystalline grains, but the need of increasing the memory densities forces the 
development of non-conventional media, like patterned media and self-assembled magnetic 
nanoparticles. The idea is to replace the randomly oriented magnetic grains with magnetic 
particles or nanodots in which it is possible to control the magnetic anisotropy orientation 
El with the advantage of reducing the noise. Whatever the magnetic support, thermal 
stability is a key ingredient in magnetic media for data storage. Increasing the storage 
density implies a reduction in the size of the magnetic grains. If the particles or grains are 
too small they lose thermal stability reaching the so called superparamagnetic limit. For 
single isolated magnetic nanoparticles, modelled as giant magnetic moments, the thermal 
stability problem has been extensively studied from the theoretical point of view flU [51 [6]|, 
however, when the system is a set of interacting particles there are still many open questions 
d [81 |9l . In particular, long-range dipolar interactions are unavoidable becoming relevant 
when increasing the packing density of the magnetic moments. In this sense, it is very 
important the understanding of the role of dipolar interactions in the relaxation dynamics of a 
set of magnetic moments due to thermal fluctuations. Patterned alloys [2] and two dimensional 
self-assembled magnetic nanoparticle arrays are systems suitable for the experimental study 
of thermal stability, since it is possible to tailor the magnetic anisotropy orientation of 
the magnetic nanoparticles with respect to the array plane. In addition, other important 
properties like the shape and size of the magnetic units can also be controlled. A landmark of 
magnetic relaxation in particle arrays is the slow, quasi-logarithmic decay of the remanent 
magnetization over several decades in time. The origin of this slow relaxation is still 
controversial. Recent simulation studies in three dimensional systems point to the necessity 
of considering very simplified models, with minimum ingredients, as a way to understand the 
mechanisms responsible for the slow relaxation and glassiness in dipolarly coupled magnetic 
particles [fTOl ITTL IT2l IT3l IT4l| . The structure of the crystal lattice (or the absence of structure), 
the nature of anisotropics (whether random or not), the presence of polydispersity in the 
particle volumes, and the volume concentration of the array, all have to be carefully considered 
in order to better understand the origin of magnetic behavior in a particular array. 

The relaxation behavior of two dimensional arrays of monodispersed magnetic moments 
has been studied theoretically [IT5l [16l [8l [9l |T7l [T8]| . Mean field calculations flTSj show that the 
relaxation may be slow without the need of introducing disorder, as commonly assumed. The 
mean field approach does not take into account correlations between particles due to dipolar 
interactions and rapidly fluctuating thermal fields [[9]|. In fact, when dipolar interactions are 
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considered, there are ordered states at low temperatures. Taking as example an array of 
particles in a square lattice and depending on the ratio between the dipolar and anisotropy 
strength, the system can order antiferromagnetically in the plane or out of the plane lfT9ll20ll . 
The complexity of the interacting problem requires the use of numerical tools in order to go 
beyond mean field results 1(171 . 

In this work we first show, by means of Monte Carlo simulations, the phase diagram of 
a system of monodispersed magnetic particles in a square lattice with perpendicular uniaxial 
anisotropy and dipolar interactions. With this information at hand, we then use a recently 
introduced Time Quantified Monte Carlo method to characterize the relaxation dynamics of 
the system. We characterize the time dependent relaxation near thermal equilibrium, which 
happens to be very well described by stretched exponential decays in the whole time window 
and in both ordered and disordered regions of the phase diagram. From the relaxation curves 
we go on to obtain the characteristic relaxation times and compare them with the predictions 
of the Brown- Arrhenius theory. In order to explore the influence of dipolar interactions on 
the effective free energy barriers that control the relaxation process, we explore a wide range 
of intensities of the dipolar interaction. This allows us to study the limits of validity of the 
Brown- Arrhenius model of relaxation and define the basic ingredients responsible for the 
observed slow dynamics. 

2. Model and Simulations 

We consider a model of a two dimensional monodispersed array of single domain magnetic 
nanoparticles arranged in a square lattice. The particles have uniaxial anisotropy which is 
oriented perpendicular to the x — y plane of the array. The classical Hamiltonian of this 
system can be written in the form: 

H = -D Si + g g 3 ^ 3 ^h. • , (1) 

i i<j T ij 

where Si = ft, * / ' p, np are three-dimensional magnetic moments of unit length. fx np is the 
value of the magnetic moment, which depends of the volume V and the magnetization M s 
of the magnetic nanoparticle, fj, np = VM S . = r^/r^ are unit vectors on the plane 
of the array. The first term is the energy contribution of the uniaxial anisotropy pointing 
in the direction of the z axis. D = K U V with K u the uniaxial anisotropy constant of the 
system under study. The second term in the Hamiltonian is the dipolar interaction between 
the magnetic particles, g = where fi is the vacuum permeability and a is the lattice 

parameter. We consider the dipolar energy in a point dipole approximation [21 J. In order to 
obtain the equilibrium phase diagram of the system we employ the Metropolis Monte Carlo 
algorithm where the spins are randomly updated in the unit sphere |[22l . For the simulations 
of the relaxation dynamics we implemented a recently introduced Time Quantified Monte 
Carlo Method (TQMC) E3ll24ll25l which allows direct comparison with experimental time 
scales. We considered the high damping limit in order to maximize the total time span of the 
simulation. The simulations were done for square lattices of linear dimensions L = 16 and 
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Figure 1. Phase diagram obtained from Monte Carlo Simulations for a system size L=32. The 
dashed line indicates the blocking temperature for the case of no n- interacting particles. 

L = 32. In order to minimize finite size effects due to the long ranged dipolar interactions, 
periodic boundary conditions were imposed by means of the Ewald method. Details for 
the implementation of the Ewald method can be found in ref. Il26l . Once the dipolar and 
anisotropy effective fields are calculated for each magnetic moment in the lattice, they are 
updated randomly using the TQMC [|25l . 

3. Results 

3.1. Phase diagram 

Figure [TJ shows the phase diagram of the system which was obtained and discussed in the 
context of the spin reorientation transition in reference 11201 . Here we plot it in a way suitable 
for discussing the dynamical behavior. All energy scales are in units of the Boltzmann 
constant, which was set to ks = 1. The dashed line in the figure corresponds to the blocking 
temperature for our simulated time window, for the case of an array of non-interacting 
particles, i.e g/D = 0. For high temperatures the particles are superparamagnetic (SP 
phase). Upon lowering T the anisotropy and dipolar terms in the energy begin to rule the 
behavior of the system, leading to two possible phases with antiferromagnetic order: for 
small temperatures and g/D < 0.8 the film orders out of plane in a checkerboard like 
configuration ll27Tl . The perpendicular staggered magnetization shows a jump at the transition 
temperature, suggesting that the transition to the perpendicular antiferromagnetic phase is first 
order. Of course, a precise determination of the order of the transition would only be possible 
considering much larger system sizes, which is beyond the scope of this work. Our results in 
the present study all refer to this region of the phase diagram. For comparison, from an array 
of 4 nm Co particles we obtain g/D = 0.67 or g/D = 0.084, if the particles are in contact or 
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are separated by a distance equal to the diameter, respectively. Here we use the Co parameters 
of Ref. [28 J. With growing dipolar coupling, the system goes through a spin reorientation 
transition to an in-plane antiferromagnet as described in EOl . 

The knowledge of the different phases in an interacting system is essential to correctly 
interpret the relaxation behavior which will be the subject of the following sections. 



3.2. Near equilibrium relaxation 

For the dynamical simulations we employed a recently introduced Time Quantified Monte 
Carlo Method (TQMC) (23l |M |25]| which allows to map Monte Carlo time steps to real 
time scales for any values of the damping parameter usually introduced in micromagnetic 
simulations. This new TQMC is very efficient and much more easy to implement and control 
than a micromagnetic simulation. The new approach is based in a controlled map between the 
stochastic Monte Carlo dynamics and the corresponding stochastic Fokker-Planck equation. 
By means of a detailed correspondence between the two approaches to the same stochastic 
dynamics, the real time scale (in units of a characteristic damping time) and MC time steps 
become related by: 

At[r K ]=a^At[MCS], (2) 
where the damping time is given by: 

tk = — • (3) 

lH k a 

In the previous relations a = D/T, R is the size of a cone where the spin is updated, 
H k = 2K/M S , M s is the saturation magnetization and a is the damping parameter. For 
details on the method and notations see references If24l 1251 |29l . In the present simulations 
0.03 < R < 0.1 in order to optimize the time span of the dynamics for the different 
temperatures. 

In order to characterize the relaxation dynamics of the system as it approaches 
equilibrium, we measured time dependent auto-correlation functions by first relaxing the 
system for 10 5 MC steps from an initially disordered configuration of the spins. Then an 
initial configuration was stored {Si(0) } and subsequent auto-correlations were defined as: 

1 N ^ 

C(t) = -J2sm-Si(t). (4) 

iv i=l 

Statistical averages of this function for many realizations of initial conditions were done. 
Auto-correlations relax similarly to the magnetization and are related to the linear response 
function by the fluctuation-dissipation theorem. 

Our reference system is an ensemble of noninteracting particles with uniaxial anisotropy. 
This corresponds to the vertical axis of figure ([I]). In figure (O we show four correlation curves 
for this case, each curve corresponding to a low temperature in the high energy barrier regime, 
D/T ^> 1. The correlation curves show a rapid decay to a rather high value where a plateau is 
developed as the temperature is lowered. This regime corresponds to the relaxation inside the 
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Figure 2. Autocorrelations for a system of noninteracting particles. Continuous lines are 
exponential fits. 



initial basin of the energy minimum near the initial condition. At long times the particles are 
able to overcome the anisotropy barriers and a final relaxation takes place where the particles 
decorrelate from their initial state. In this case of noninteracting particles the relaxations are 
exponential in time, as can be seen in the fits in figure ©: 



Note that the horizontal time scale is given in units of (7 Hk) -1 . Taking as reference cobalt 
with a high damping constant a = 1, then tk = 4.82 x 1CT 11 s using 7 = 0.2212 10 6 m/As^ 1 
for the gyromagnetic ratio and H k = 168 k Aj 'm for the anisotropy field ||28i 

From the results for the correlation functions we can obtain the relaxation times for the 
particles. In a classic paper, Brown obtained approximate expressions for the relaxation time 
of single domain magnetic particles. In the case of zero applied field and high energy barriers 
the result of Brown is [4] : 



In figure © we show the relaxation times obtained from the exponential fits of figure © 
together with the analytic approximation by Brown in a log-linear plot. The agreement is very 
good except for the point corresponding to the lowest temperature. This departure may be due 
to insufficient relaxation of the particles for this very high energy barrier. This explanation can 
be supported by calculating the blocking temperature for the ensemble of particles, defined 
as the temperature for which the experimental time scale equals the relaxation time. This 
Tb is usually calculated assuming an exponential or Arrhenius behavior of the relaxation 
time. This is not completely true in this case, where the prefactor of the exponential in © 
depends on temperature. Nevertheless a fit with a purely Arrhenius behaviour of the form 



C(t) = C e-^. 



(5) 




(6) 
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Figure 3. Relaxation times for a system of noninteracting particles. The solid line is the 
analytic approximation by Brown. 



t = r exp (AE/T) works reasonably well and allows a rough determination of the blocking 
temperature: 

DAE 

T B = 7T , r , (7) 
In (texp/To) 

where in our simulations t exp ~ 10 6 . From the fitting values of AE = 0.85 and r = 2.27, we 
obtain T B / D « 0.065 which is near the value of the lowest temperature simulated. Then it is 
probable that the curve corresponding to T = 0.066 in Fig. © be still subject to considerable 
fluctuations. 

Upon switching the dipolar interaction the main observation is that the relaxation is no 
more exponential. In figure @] we show correlation curves for a small intensity of the dipolar 
interactions g/D = 0.0625 and different temperatures, all in the superparamagnetic phase 
(see figure [[])• Together with the data points, the solid lines represent fits with a stretched 
exponential form: 

C{t) = C e~W. (8) 

In this case of weak dipolar interaction the values of the exponent (3 are around 0.9. In fact 
we found that the typical (3 values diminish as the intensity of the dipolar coupling grows. 
Moreover, the stretched exponential form can describe the whole time span of the relaxation 
for any value of the dipolar intensity. Nevertheless, to our knowledge, there is no analytic 
prediction for this particular form of the time dependent relaxation in systems of interacting 
particles. 

A quasi-logarithmic time decay is widely employed to fit relaxation data in systems with 
some kind of disorder. Nevertheless some time ago Dahlberg et al. [1301 showed that this slow 
relaxation need not be associated with disorder in the system and can as well be present due 
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Figure 4. Correlations for g/D — 0.0625. Solid lines are fits to stretched exponential forms 
(see text). 

to the long ranged dipolar interactions. They verified this by solving numerically a mean field 
model of two state particles. More recently Denisov et al. JH obtained analytic expressions 
for the relaxation of the remanent magnetization in a 2d system of interacting dipoles with 
strong perpendicular anisotropy from an initially saturated state. They find a crossover from 
slow relaxation at intermediate times to exponential one at very long times. Although the 
models are similar, a direct comparison of their results with ours is not straightforward due to 
the different protocols used in both cases. 

The stretched exponential relaxation is observed also in the low temperature regime of 
the present system, when it relaxes towards a state with long range antiferromagnetic order 
out of the plane. In figure © we show autocorrelation curves for temperatures corresponding 
to the ordered phase. For this value of g/D = 0.25 the stretching exponent (3 is around 0.6, 
i.e. the relaxation is much slower than in the non-interacting case. Nevertheless, although 
not shown in the figure, the behaviour of relaxation times does not present any anomaly when 
crossing the phase transition temperature. Due to the first order nature of the transition (and 
the small size of the system) only a gradual change of stability between the disordered and 
ordered phases takes place, with no evident effect on the typical relaxation times. 

Together with the departure from exponential decay, when interactions are switched 
on the relaxation times also behave differently than the independent particle relation eq. 
©. In figure © the characteristic relaxation times corresponding to dipolar to anisotropy 
ratios g/D = 0.0625,0.125 and 0.25 are shown. The values of r were obtained from the 
corresponding stretched exponential fits and the data was fitted with the functional form: 

— K^(aAE)- 1 / 2 e a * E , (9) 
t k 2 
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Figure 5. Correlations for g/D = 0.25. Solid lines are fits to stretched exponential forms (see 
text). 
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Figure 6. Relaxation times for a system with three different dipolar to anisotropy interaction 
ratios. The solid lines are fits with expression (0 as described in the text. 
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which is a natural modification of Brown's relation ©, where now DAE represents an 
effective free energy barrier (we fixed D = 1 in the simulations whenever g ^ 0). Typical 
values of AE are 0.99, 0.83 and 0.45 for g/D = 0.0625,0.125 and 0.25 respectively. It 
was not possible to fit the data for the interacting systems fixing AE to one, as in the 
noninteracting case, even considering the exponent of the temperature dependent prefactor 
as a fitting parameter different from 1/2. It is important to note that while the data for 
g/D = 0.0625 all correspond to the superparamagnetic phase, those corresponding to the 
two larger values of dipolar interaction probe the system relaxing towards an ordered phase. 
In spite of this, it can be seen in figure © that in all cases the data is compatible with 
a single free energy barrier as represented by the nearly Arrhenius form ©. From this 
behavior it can be concluded that glassy behaviour is not present in the dipolar interacting 
system, as this case should imply a non-Arrhenius dependence of relaxation times. Of course, 
some kind of disorder not considered in the present simulations, e.g. random anisotropics or 
polydispersity, may lead to spin glass like behaviour (see, for example, discussions in |[TTl[T2l 
and references therein). Instead, we observe that dipolar interactions are responsible for an 
effective reduction of the free energy barriers which control the relaxation process, even with 
an initial state where demagnetization fields are significantly reduced. This energy barrier 
reduction is related to the growing of the ordered state as can be realized by looking at the 
order parameter evolution during the relaxation. It is worth to stress that the growing of the 
ordered phase does not affect the Arrhenius character of the relaxation process. 

4. Conclusions 

We have studied the relaxation dynamics of an array of dipolar interacting magnetic particles. 
Instead of looking at the remanent magnetization from a saturated state, we have first let 
the system relax near equilibrium and the subsequent relaxation was studied. Our model 
of a monodispersed array of single domain particles with uniform perpendicular anisotropy, 
while clearly simplified, is a first step towards a more systematic study of the origin of slow 
dynamics in realistic nanoparticles arrays. We obtained two main results: first, regarding the 
time dependence of relaxation functions, we found that while the relaxation is exponential 
in the noninteracting case, it changes to a slower relaxation form as the dipolar interaction is 
switched on, already for a very weak interaction strength. A stretched exponential form fits 
remarkably well the data over the whole time span of the relaxation, contrary to the known 
limitations of a logarithmic fit Ifi~5ll31~ll32~1l . While this may indicate a fundamental role of the 
stretched exponential form as a consequence of the microscopic interactions [fT0l[TTl[T2l . this 
remains to be proved on theoretical basis. We also found that the slow dynamics is present 
in the interacting system regardless of the temperature regime, while relaxing towards the 
superparamagnetic or ordered phases. Second, we showed that the Brown- Arrhenius result for 
the temperature dependence of relaxation times works very well for independent particles, but 
needs to be modified when interactions are present. We proposed a slight modification of the 
original Brown's results in order to fit the data. Within this form, the behaviour of relaxation 
times shows that effective free energy barriers are reduced with growing dipolar interactions, 
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at least in the present case where the system is initially already in a demagnetized state. 

The Time Quantified Monte Carlo method turned out to be a convenient and efficient way 
to perform simulations which can be directly confronted with experimental data. We used the 
example of cobalt nanoparticles in which we showed that realistic time scales can be reached 
within reasonable model parameters. 

Monte Carlo simulations of three dimensional arrays of particles, including both 
configurational and anisotropy disorder IfTOl ITT1 [T2l [P3l [141 . have shown that anisotropy 
disorder smears out the effect of configurational disorder |fT2|. and that this disorder 
is particularly important when the ground state is a columnar antiferromagnet ffT2l 
[T3l . Moreover, glassy behavior, present in these cases, seems to be enhanced in the 
configurationally ordered case lfl4ll . Then, future studies of two dimensional systems should 
add the effects of dipolar interactions with the influence of disorder on relaxation, present 
for example when the anisotropy axes are randomly oriented, which is the usual case in 
self-assembled or epitaxial nanoparticle systems lT3~2l [331 l34ll . For the important case 
of perpendicular anisotropy, relevant for data storage, disorder can also be present in the 
intensities of local anisotropics ll35l l34l . Polydispersed assemblies with a distribution of 
grain volumes can also be straightforwardly considered within the present model ll32l [33l . 
We plan to address this issues in future works, together with different dynamical protocols 
more directly related with experiments, as for example, the behaviour of the thermoremanent 
magnetization in field cooled protocols and the influence of external fields after a zero field 
cooled process. 
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